Melting of trapped few particle systems 
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In small confined systems predictions for the melting point strongly depend on the choice of 
quantity and on the way it is computed, even yielding divergent and ambiguous results. We present 
a very simple quantity which allows to control these problems - the variance of the block averaged 
interparticle distance fluctuations. 
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Crystallization and melting and, more generally, phase 
transitions are well known to pertain to very large sys- 
tems only. At the same time, solid-like or liquid-like 
behavior has been observed in finite systems containing 
only one hundred or even 10 particles and is becoming of 
increasing interest in many fields of physics, chemistry, 
and beyond. Current examples include bosonic crystals 
and supersolids, e.g. electrons or excitons in quan- 
tum dots 0, ions in traps [1], dusty plasma crystals [1], 
atomic clusters [1, @|, polymers 0] etc. The notion of liq- 
uid and solid "phases" has been used successfully to char- 
acterize qualitatively different behaviors which resemble 
the corresponding properties in macroscopic systems and 
will be used here as well, following the definition of ref. 
. From the existence of phase-like states in very small 
systems arises the fundamental question of how to char- 
acterize phase changes and further, how many particles 
does it take at least to observe a phase transition. 

In macroscopic systems a solid-liquid transition can 
be verified by a variety of quantities including free en- 
ergy differences, order parameters, specific heat, trans- 
port properties, structure factors, correlation functions 
and so on, e.g. which yield more or less equivalent 

results for the melting point. A particularly simple and 
transparent quantity is magnitude of the particle posi- 
tion fluctuations normalized to the interparticle distance 
(Lindemann ratio ml), for an overview see. But 
when applied to two-dimensional (2D) systems, ul shows 
a logarithmic divergence with system size This led 
to modified definitions, including the relative interparti- 
cle distance fluctuations (IDF) [13, [3, [T^ 
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which are also well behaved in macroscopic 2D and ID 
systems. Here = jr^ — rj| is the distance between 
two particles and (. . .) denotes thermal averaging. In 
macrosopic systems u^d shows a jump at the melting 
point which clearly reflects the increased delocalization 
of particles in the liquid phase compared to a crystal. 

However, when applied to small systems, N < 100, nei- 
ther Ml nor Mrei exhibit a jump upon classical or quantum 
melting, but rather a continuous increase over some finite 



temperature or density interval ~ a familiar finite size 
effect. Therefore, it is very difficult to determine a tran- 
sition point and the critical magnitude of the fluctuations 
M^g}*. Even worse, the result for Mroi (and hence the melt- 
ing point) depend crucially on the method of calculation 
and on its duration. Increasing the length of a simulation 
(and the expected accuracy) may lead to growing system- 
atic errors predicting a too low melting temperature, as 
was noted by Frantz Q and a few others. This is, of 
course, critical for reliable computer simulation of phase 
transitions in finite systems. In this Letter, we analyze 
the reasons of this behavior and present a solution. We 
propose a novel quantity, the variance of the block aver- 
aged interparticle distance fluctuations, which is sensitive 
to melting transitions and does not exhibit the conver- 
gence problems of Mrei- We demonstrate the behavior of 
this quantity both, for classical and quantum melting by 
performing classical Monte Carlo (MC) and path integral 
Monte Carlo (PIMC) simulations, respectively. 

Model and parameters. While our approach is gen- 
erally applicable we concentrate on strongly correlated 
classical or quantum particles in a parabolic trap in 2D 
and 3D described by the Hamiltonian 
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The system is in a heat bath with temperature T and 
has a fixed particle number N (canonical ensemble). Be- 
low we use the dimensionless temperature k^Tro/e^ —> T 
where rg denotes the ground state distance of two par- 
ticles, Tg = 2e^ jmuP'. For quantum systems, the cou- 
pling parameter is A = e^ j^l^hiS) where Zq is the oscillator 
length Zg = fij (moj). The ground state of this system con- 
sists of concentric spherical rings (2D), cf. Fig.[T]or shells 
(3D) d, Q. This model has been very successul in de- 
scribing trapped particles in many fields and has the ad- 
vantage that classical melting (by temperature increase) 
and quantum melting (via compression by increasing lu), 
including spin effects [3, [3], can be analyzed on equal 
footing |16| . 

Liquid and solid "phases". The potential energy 
landscape of the system ([2]) has numerous local minima 
but, in contrast to other finite systems such as atomic 
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FIG. 1: Configuration of a 2D trapped quantum system of 
N = 8 spin-polarized bosons described by Eq. ([2}. Fig. a) 
shows the liquid state (Ai — 14) and b) the solid state (A4 — 
30), temperature is close to the ground state. Figure c) shows 
the radial density profile p(r) for both configurations. 

clusters 0, 01 , they are not associated with "phases" but 
rather correspond to the ground state and metastable 
states (e.g. different shell configurations) which often are 
energetically very close, e.g. [Tj). With increasing tem- 
perature an increasing number of these states becomes 
occupied. Melting proceeds as an isomerization transi- 
tion with the system switching rapidly between a fast 
growing number of different configurations above some 
threshold temperature p^ . 

However, with reduction of N the number of station- 
ary states decreases until only the ground state remains. 
This is the case for iV = 4 in 2D which, due to its simplic- 
ity, allows for a transparent analysis of melting processes 
m the system The pair distances show a character- 
istic behavior as a function of simulation time (MC step) 
k, cf. Fig. [^a): oscillations around some average value 
followed by a jump to a different value and again oscilla- 
tions around a different mean and so on. This is readily 
understood: in its ground state the particles occupy the 
corners of a square of length a, so there exist two possi- 
ble values for the six pair distances: a and the diagonal 
b = ^/2a, which are the mean values around which the 
distances fluctuate. A jump occurs whenever two parti- 
cles i and j exchange their positions. Then the distances 
rik and rjk to the remaining particles will change. While 
this leads to the same ground state (permutational iso- 
mer) this process costs energy associated with overcom- 
ing of a potential energy barrier. With increasing tem- 
perature, the frequency Vj of these jumps grows steadily 
until around T = T2 a rapid growth of Vj is observed. 
Finally, aX T = T^, pair exchanges occur constantly (bot- 
tom of Fig.[2]a), and particles are practically delocalized. 
This behavior of Vj clearly resembles a "phase transition" 
with the melting point being located inbetween the two 
limits Ti (solid) and T4 (liquid). 

We verify this hypothesis by computing the IDF, 
Eq. ([1]) for this system, cf. Fig. [Hb). At low temper- 
atures, Urol is small, slowly increasing with T. Around 
T = T2 the increase steepens slightly (rightmost curve). 
Repeating the calculations with higher accuracy, by sub- 
sequently increasing the simulation length L (number of 
MC-steps) by factors 10, 100, 1000, Urci shifts left towards 
smaller temperatures, and no convergence is observed. 
Thus, longer calculations yield an increase of u^d already 




MC steps 



FIG. 2: (a) Distance of an arbitrary pair of A'" = 4 classical 
particles in 2D as a function of MC step. From top to bottom: 
Ti = 0.02 (solid-like), T2 = 0.06 and T3 = 0.09 (transition re- 
gion), and Ti = 0.5 (liquid-like), a and h = \/2a denote 
the two possible interparticle distances in the ground state, 
(b) Temperature dependence of the mean block averaged IDF 
Urcu for different block lengths M = 10^ lO", 10^ 10*^ (right 
to left) [equivalent to computing Uroi, Eq. ([l|, from multiple 
simulations of length L — M]. (c) The corresponding second 
moment cr^^^j, Eq. ((3}. (d) Specific heat Cv and energy cor- 
relation time fccorr. (e) Total energy autocorrelation function 
Ce, Eq. Q, for three temperatures. 



in the solid-like regime, even though jumps are very rare, 
so the results for Uj-ci, Eq. ([T]), are ambiguous and un- 
reliable. The reason is that, even in the solid state, a 
jump will be captured if L is sufficiently long. This im- 
mediately leads to a significant increase of Uioi emulating 
liquid- like behavior d9j . Similar observations were made 
for clusters in Ref. @. 

Solution of the convergence problem of Urei- We 

solve this problem by sub-dividing the time sequence in 
K blocks of equal length M {L — K-M) and compute the 
block averaged IDF u-cc\{s) according to Eq. ([1]) for each 
block s [13, [2l| and its mean Uici = SfLi '^ncii.s). 
To suppress the influence of jumps to Urci in the solid 
regime, M must be chosen small enough to restrict jump- 
related contributions to a small number of blocks and, 
at the same time, large enough to allow for convergence 
of contributions related to local vibrations. This choice 
does not influence the convergence of Urei in the liquid 
regime which is dominated by frequent jumps on a time 
scale comparable to that of local vibrations and, hence. 
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FIG. 3: Left: Typical behavior of the block averaged IDF vs. 
block number s for N = 8 charged bosons in 2D for different 
coupling strengths A: Ai = 14, A2 = 22, A3 = 26 and A4 — 30. 
Each point is an average over a block of length M — 1000. 
Right: Histograms show the probability P of different values 
Mroi averaged over a total of 9000 blocks. Results are from 
PIMC simulations with 200 .. . 500 beads HI] of system 



well below M. We demonstrate the behavior of Uici(s) 
for a quantum phase transition of iV = 8 bosons in 2D, 
cf. Fig. [3l In the solid regime there are rare spikes of 
itrci(s) corresponding to occasional blocks containing one 
jump leading to a sharply peaked probability distribution 
P(urci). In the transition region, however, each block 
may "catch" from zero to a few jumps, so the fluctuations 
of Urci(s) increase and P(uici) broadens. Finally, in the 
liquid regime, jumps occur with an almost constant rate 
in every block, so the fluctuations of Urei(s) are small 
[P(mi-ci) has again a single sharp peak], while the mean 
is shifted to a higher value above 0.3, typical for a liquid. 

From this we conclude that, in the vicinity of the 
melting transition, the width of the distribution P{urci) 
reaches a maximum. This behavior is well captured by 
the second moment of Mioi(s), i.e. the variance of the block 
averaged interparticle distance fluctuations (VIDF) 
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This allows us to obtain a reasonable estimate of the 
melting temperature T^"^'* from the peak of o^^^^ (T) [l^l ■ 
Note that Mrei is sensitive to the jump frequency v^^ in 
contrast to itrei of Eq. ([T]). The sensitivity does depend 
on the block length M : larger M cause an increase of 
Wrci (as discussed before) and shift the maximum of cr„^^| 
to lower temperatures, cf. Fig. [2lb)-c). 

Therefore, to properly choose M an independent quan- 
tity is needed which should not require block averaging 
and be invariant with respect to particle exchanges and 
pair distance jumps. A quantity fulfilling these require- 
ments is the total energy E and its autocorrelation func- 
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FIG. 4: Mean value Ujoi (top) and second moment a^.^^^ (bot- 
tom) of the block averaged IDF for different particle num- 
bers A'^. Left: temperature dependence of a classical 3D sys- 
tem (classical melting, classical MC simulations). Right: 2D 
quantum system, dependence on the quantum coupling pa- 
rameter A (quantum melting, PIMC results). In both cases 
the block length equals M = 1000. Dashed lines locate the 
critical values of T (or A) and tire!'- 
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We found that the decay rate of Csik) varies non- 
monotonically with temperature where the slowest decay 
is observed just in the transition region, cf. the example 
shown in Fig. [He). This suggests that the correlation 
time, fccorr(r) = EfcC'ij(fc,T), cf. Fig. [Hd), is sensitive 
to thermal melting, allowing us to identify the melting 
temperature TJf'' from the maximum of fccorr- Compar- 
ing the values T^^^ and r"''*(M) provides a straight- 
forward way to identify the proper block length M . In 
all cases of thermal melting we investigated agreement 
is found for M in the range of 1000 . . . 10000, where the 
common definition of a Monte-Carlo step is used [23j . 

We mention that in the case of quantum melting the 
situation is more complex. Nevertheless, we found that 
the same range of M seems appropriate here as well, 
however, the analysis requires to use a combination of 
different quantities such as the pair distribution or bond 
angular symmetry parameters etc. 

Applications. We have verified the behavior of the 
VIDF, (7tij.^j, for a large variety of classical and quantum 
systems described by Eq. ^ of various sizes and dimen- 
sionality. As a first illustration we show in Fig. 2] (left 
side) MC results for a classical 3D system of = 4... 20 
particles the state of which is completely characterized 
by the temperature T. One clearly sees that in all cases 
"Urol increases with T, but for small A^ the reduction is 



4 



very gradual, not allowing us to single out a "melting 
temperature" from u^ei- At the same time, in all but 
one case cru,oi has a well pronounced peak at a certain 
T which is identified as T'^'''*. Also, the critical value of 
the fluctations may be deduced from the peak position of 
cr„,,^j yielding uj:^}* « 0.08 ... 0.16 which is in good agree- 
ment with macroscopic classical Coulomb systems. Note 
the special case of iV = 5 showing a low value of T'^'''* 
which is well known and explained by the low symmetry 
of this cluster [l3] ■ While this behavior is hardly visible 
in Uj-ei it is clearly detected by cr„^^j . 

As a second example we consider quantum melting 
upon compression in a 2D system of spin polarized 
charged bosons at very low temperature close to the 
ground state. Calculations for particle numbers up to 
iV = 60 were done using PIMC simulations, for details 
see e.g. [l^. Right hand side of Fig. [4] shows results for 
three cases, N = 19,20,21, more examples are given in 
Ref. For large A, the particles are localized resem- 

bling a crystal as seen in Fig. [TJ Decrease of A is asso- 
ciated with increasing wave function overlap and eventu- 
ally quantum melting by tunneling of particles between 
lattice sites. Again we observe a gradual reduction of MicI 
when A is increased. In contrast, Uu^^j has a pronouced 
peak which allows us to determine the critical value of A 
to A ~ 25 ... 30 depending on the particle number. The 
corresponding critical fluctuations, « 0.22... 0.25, 
are again close to the value known from simulations of 
macroscopic Bose systems. 

These two examples are representative for the classical 
and quantum melting behavior of the system ^ , also for 
other pair potentials. All our calculations have confirmed 
the robustness and efhciency of the VIDF for the anal- 
ysis of melting in small systems. We can now proceed 
and analyze the question what is the minimum system 



size to observe crystallization or melting? Our simula- 
tions have revealed that du^^i has a maximum for parti- 
cle numbers as small as 4 in 2D and 5 in 3D. In contrast, 
for 4 particles in 3D, cr„_,^,| shows a monotonic increase, 
see Fig. |4] (top left). This is easily understood. The 
ground state of 4 (3) particles in 3D (2D) resembles an 
unilateral tetraeder (triangle) and has only a single in- 
terparticle distance. Thus, a jump (pair exchange) does 
not alter the distribution of pair distances, and cru,ei ti^s 
no maximum. 

In summary, we have proposed a novel quantity - the 
variance of the block averaged interparticle distance fluc- 
tuations - which is sensitive to fluctuations in finite sys- 
tems. A maximum of (Ju,^^i allows one to reliably detect 
the existence of structural changes which are analogous 
to solid-liquid phase transitions in macroscopic systems. 
It further directly yields a consistent estimate of the melt- 
ing point [22*] and the critical fluctuations u^l*" in classi- 
cal and quantum systems, thereby curing the sensitivity 
and convergence problems of the conventional distance 
fluctuation parameters. While for classical systems the 
energy autocorrelation function Ce allows for a calibra- 
tion of the block length, this does not work for quantum 
melting where further analysis is required. Also, it re- 
mains an interesting question to analyze the behavior of 
(7„^,^j in other finite systems, including atomic clusters or 
homopolymers etc, as well as in time-depcndend simula- 
tions (such as molecular dynamics). Finally, in the case 
of strongly inhomogeneous macroscopic systems where 
melting is known to proceed via a sequence of different 
processes, the VIDF should allow for a deeper insight and 
a space-resolved analysis of the fluctuations. 
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